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Abstract 



Motivation: 



Many biochemical pathways are known, but the numerous parameters required to correctly 
explore the dynamics of the pathways are not known. For this reason, algorithms that 
can make inferences by looking at the topology of a network are desirable. In this work, 
we are particular interested in the question of whether a given pathway can potentially 
harbor multiple stable steady states. In other words, the challenge is to find the set 
of parameters such that the dynamical system defined by a set of ordinary differential 
equations will contain multiple stable steady states. Being able to find parameters that 
cause a network to be bistable may also be benefitial for engineering synthetic bistable 
systems where the engineer needs to know a working set of parameters. 



Result: 



We have developed an algorithm that optimizes the parameters of a dynamical system so 
that the system will contain at least one saddle or unstable point. The algorithm then 
looks at trajectories around this saddle or unstable point to see whether the different 
trajectories converge to different stable points. The algorithm returns the parameters 
that causes the system to exhibit multiple stable points. Since this is an optimization 
algorithm, it is not quaranteed to find a solution. Repeated runs are often required to 
find a solution for systems where only a narrow set of parameters exhibit bistability. 



Availability: 



The C code for the algorithm is available at http:/ /tinkercell. googlecode.com 



1 Introduction 



The existence of more than one stable state is a common theme in biochemical networks. 
Bistable networks, or networks with two stable steady states, have been observed to be 
very common in biology and may have great importance in the overall behavior of a 
single cell or multi-cell system. Theoretical methods have been developed to determine 
the possibility of bistability by looking at network topology [TJ El HI El E] . However, 
these methods are usually limited to networks that obey mass-action kinetics, but a vast 
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majority of models use more complex kinetics. Other theoretical methods also exist that 
use certain criteria to guarantee the presence of multiple stable states |7j. These meth- 
ods are based on the understanding that a positive feedback is a neccessary condition for 
bistability [8J |9j, and monotonicity, in addition to positive feedback, can be a sufficient 
condition. However, such methods can be stringent, thus missing some of the bistable sys- 
tems. A numerical method has also been developed that uses the eigenvalues of a system 
to search for bifurcation points, thus identifying the point of transition from monostability 
to bistability [ID]. The algorithm presented in this work is also a numerical method, but 
it uses a slightly different method to optimize the system which shows improvement over 
the aforementioed algorithm. 



2 Method 



The goal of the algorithm presented in this work is as follows: 



Given: A system of differential equations, dX/dt = F(X,P) with an unknown 
parameter set P = ki...k n 

Find: P such that dX/dt = F(X, P) has more than one stable steady state. We define 
a stable steady state as the point where a time-course simulation converges and where 
all the eigenvalues of the Jacobian matrix are negative. 



We postulate that a system with more than one stable point will have at least one unstable 
or saddle point. The argument for this statement can be understood if one considers the 
energy landscape of a dynamical system. The energy landscape is a hypothetical phase 
plane where the "energy" at each point is proportional to the probability of the system 
existing at that point at any given time. The valleys, or low energy regions, of such an 
energy landscape correspond to the stable points, and the peaks, or high energy regions, 
correspond to the unstable or saddle points. One can argue that if there are two valleys, 
there will exist a peak between them. Therefore, a system with two or more stable points 
will have an unstable or saddle point between them. However, one should note that the 
existence of a peak does not acertain existence of two valleys. However, for the algorithm 
presented in this work, we assume that the majority of the cases for biological systems 
will be such that the existence of a peak also implies the existence of two or more valleys. 
Under this assumption, one can find bistable systems by finding systems with one unstable 
or saddle point. This is the goal of the algorithm presented in this work: search for a 
system that contains an unstable of saddle point. Once this point, or peak, is found, 
the two stable points, or valleys, can be found by following the trajectories around the 
unstable or saddle point. 
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Given the above general idea, one can restate the goal of the algorithm differently: find 
P such that dX/dt = F(X,P) has at least one unstable or saddle point. 

Once the unstable or saddle point is found, the existence of more than one stable point 
can be confirmed through time-course simulations around the unstable or saddle point. 
In order to find the parameters that cause the system to harbor at least one saddle or 
unstable point, the algorithm optimizes the following objective function: 



Given: A system of differential equations, dX/dt = aF(X, P) with an unknown 
parameter set P = k\...k n and vector a. a is a vector of reals and has the same size as X. 

Find: P and a such that not all values in a are positive and dX/dt = aF(X,P) has at 
least one stable steady state. 



Since the a vector described above is not all positive, the stable steady state of the 
system dX/dt = aF(X,P) will be a saddle or unstable steady state for the original 
system dX/dt = F(X,P). Note that the steady state solutions for the system with or 
without a is the same since a is a vector of constants. However, since at least one of 
the a values is negative, the stability of each steady state will be different for the system 
with and without a. For a monostable system, the system with a will not converge at 
all, because the only steady state in the system will be a saddle or unstable point. For 
a bistable system, there will be some a such that the system with a will converge to the 
saddle or unstable point of the original system. Therefore the parameter set, P, that 
causes dX/dt = aF(X,P) to converge will cause dX/dt = F(X,P) to be bistable. 



2.1 Numerical libraries used 



We used a genetic algorithm to optimize the objective function. The genetic algorithm 

C library was obtained from |http: / /tinkercell.googlecode.com} The CVODE library in 

the Sundials package was used to perform numerical integration. The optimization was 

performed using a Nelder-Mead C library obtained from 

http: / / www.ritsumei.acj p/se/ ~hirai / edu / 2003 / algorithm/index-e.html 

The eigenvalue calculations were performed using CLAPACK. 
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3 Results 



We tested the algorithm using a few known bistable systems. One of the systems that 
was tested was: 



02(^4 



k^x + kox 2 y — k\x 3 ) 
fay ~ (kox 2 y - hx 3 )) 



The optimization returned the following values: where the a 
k = [4.27, 2.39, 0.007, 5.55, 4.44, 0.5] 



-0.3995,0.917] and 



The system above is shown in Figure [Ta] and [TbJ The a and k stated above were found by 
forcing the system to have one stable state. Note that a has at least one negative value. 



The stable point in Figure la is the saddle point in lb which will always be the case. 
The phase portrait in Figure lb is the same system as Figure la except with a = [1,1], 
i.e. the normal system. The normal system with the parameter set k is bistable because 
the system with a has a stable point. 



Another system that was tested is a genetic toggle switch [TT] : 



x = ai(k0/ (ki + y 4 ) - k 2 x) 
y = a 2 (k3/(k 4: + x 4 ) - k 5 y) 



The optimization returned the following values: where the a 
[2.8,0.5,1.04,2.2,0.25,1.05] 



•1.0,2.5] and k 



Again, since the system with a has a stable point, the system without a has a saddle 
point and two stable points. The phase portraits are shown in Figure [2} 



The algorithm is not always successful in finding the correct parameters. For the systems 
described above, the algorithm is generally successful once in every two runs. It takes 
roughly 20 iterations to find the parameters, which is an average of 20 seconds (1 second 
per iteration) on a computer using AMD Opteron 165 1.8GHz. 



Previous methods were unsuccessful to identify bistability in both of the systems. The 
second system (Figure [2]) does not use mass-action kinetics, so methods that rely on 
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(a) unstable steady state 



x ' = 1 (k2 - k3 x + kO x x y - k1 x x x) kO = 4.27 k1 = 2.3S k2 = 0.007 

y 1 = 1 (k4 - k5 y - (kO x x y - k1 x x x)] k3 = 5.55 k4 = 4.44 k5 = 0.5 



D 02 04 0.6 08 1 1 2 1.4 1 6 1 8 2 



(b) stable steady states 

Figure 1: (a) The phase portrait of a 2-variable system when one of the a values is 
negative, (b) The phase portrait of the same variable when all of the a values are 1. Note 
that the stable point in (a) is the unstable point in (b). 
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x ' = - 1 {kO/(k1 + v 4 ) - k2 x) 
y ' = 2.5 (k3/{k4 + x 1 ) - k5 y) 




(a) unstable steady state 



x ' = 1 (kO/(k1 + y 4 ) - k2 x) 
y ' = 1 (k3/(k4 + x 1 ) - k5 y) 




(b) stable steady state 

Figure 2: (a) The phase portrait of a 2-variable system when one of the a values is 
negative, (b) The phase portrait of the same variable when all of the a values are 1. Note 
that the stable point in (a) is the unstable point in (b). 
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mass-action kinetics will fail. The first system (Figure [TJ is not always monotonic, so 
using monotonicity and positive feedback as a predictor for bistability will fail. Previous 
numeric methods (ref) were also unsuccessful in finding correct parameters for either of 
the systems. 



4 Conclusions 



We have presented an algorithm that can find parameters for a system such that the system 
will contain more than one stable steady state. The application for such an algorithm 
would be for systems, particularly biological systems, where parameters are generally not 
known. Using this algorithm, one can "guess" that the system can potentially be bistable. 
If the bistability of the system holds biological meaning, then it is possible that the system 
is indeed bistable. 



The algorithm presented in this work is not guaranteed to find the parameters for bista- 
bility, which can be an issue for systems where the range of parameters for bistability is 
small. However, we hope that the general concept from this algorithm can be modified 
and improved. For example, the genetic algorithm that is used for optimization can be 
replaced with a different algorithm. The objective function can be restated to improve 
the search process. Such changes may improve the success rate of the algorithm. 
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